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ABSTRACT 


The following is a survey article designed to provide an intro- 
duction to the subject of turbulence modeling, and to explain the need 
for such models. 

The subject is developed along chronological lines since this 
provides a logical development plan and also because it then moves 
from relatively simple phenomenological models through more compli- 
cated procedures and ultimately to the subject of large-eddy simulation. 



"...a numerical procedure without a turbulence model 
stands in the same relation to a complete calcula- 
tion scheme as an ox does to a bull." 

Peter Bradshaw 

I. INTRODUCTION 

For over a hundred years intelligent people have worked long hours 
and written thousands of technical papers on the last unsolved problem 
of classical physics, turbulent fluid flow. Indeed, there appears to be 
an emerging belief that the problem may not have a solution in the usual 
sense of the word. That is to say that turbulence seems to be simply 
the manifestation of ensemble or time averages of bounded instabilities 
resulting from multiple bifurcations in solutions of the Navier-Stokes 
equations. This is not exactly a helpful phenomenological explanation, 
but it may well be the best we will be able to give. 

Because of the absence of a consistent and comprehensive theory of 
turbulent flows, their prediction cannot be made from first principles, 
but must be based on semiempirical models. Modeling the physics rather 
than solving the full, unsteady, three-dimensional Navier-Stokes equa- 
tions is imposed by the fact that even the largest of today’s computers 
are inadequate to predict any real flow. 

In this chapter, rather than attempting to treat the multiple 
potential applications of turbulent-flow predictions, we will consider 
the background and development of our current predictive capability in 
a more or less general setting. Several other works presenting more 
specific information, or an alternative point of view, can be found in 
references [ 1 ] - [5 ] . 
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Hinze [1] defines turbulence as follows: "Turbulent fluid motion 

is an irregular condition of flow in which the various quantities show 
a random variation with time and space coordinates, so that statistically 
distinct average values can be observed." A few minutes of observing 
smoke emerging from a tall stack will be sufficient to provide a feeling 
for the above definition. 

There is a consensus among researchers of turbulence that the 
unsteady, three-dimensional Navier-Stokes equations solved on a fine 
enough mesh can describe turbulent flow from first principles. Eventu- 
ally, given computers which are sufficiently large, sufficiently fast, 
ar.d sufficiently cheap, we will be able to solve these equations 
directly, and the problem of turbulence modeling will simply vanish. 

Thus, the need for turbulence modeling is caused solely by our inability 
to solve the equations that describe the physical problem on a fine 
enough mesh to resolve all the relevant physical scales. The obvious 
question which arises is, why must anyone continue to work on it if, in 
t ‘.me , the problem will go away? The answer is that the time scale for 
tae development of the requisite computational power, on machines and 
algorithms, is very large. Lomax [6] has pointed out that existing 
computer memories will accommodate about a 64 mesh for compressible 
turbulent-flow calculations, whereas a computer memory on the order of 
512 3 mesh (roughly a three-order-of -magnitude increase) would be required 
to compute a portion of an incompressible, homogeneous, turbulent flow. 
Furthermore, for any industrially significant flow, meaningful calcula- 
tions of turbulence from first principles would probably require at 
least another two-order-of-magnitude increase of computer memory. Such 
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a powerful computer would permit computation of eddies of a scale larger 
than that associated with the Kolmogorov equilibrium scale, coupled with 
a hypothesized universal model at smaller scales. Since the required 
machine time varies typically with the square of the number of mesh 
points, most likely 10 orders of magnitude increase in computer power 
is needed before we are in a position to handle real turbulent-flow 
problems from first principles on a routine basis. Without a break- 
through in either computer design or turbulence theory it would appear 
that the millennium is 50 to 100 years in the future. This bleak, 
long-term outlook provides the motivation to seek turbulence models for 
application to the Reynolds -averaged equations. These equations will be 
discussed in detail in subsequent sections. 

The techniques of turbulent-flow prediction currently in use or 
under development follow traditional lines since much of the foundations 
of current theory were developed long before the digital computer came 
into significant use. To understand these traditional approaches, let 
us return to Hinze’s definition of turbulence. The relevant elements, 
for our purposes, can be restated as: the flow is irregular, unsteady, 

three dimensional, and statistically distinct averages exist. It is 
this last property which makes the flow approachable at all. 

In the following sections we will describe the equations of turbu- 
lent flow and show how the need for modeling arises, discuss some of the 
representative examples of today’s turbulence models and some aspects 
of the computational schemes required, touch on some special topics 
(i.e., separation, unsteady flows, and three-dimensional flows), and. 
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finally, briefly consider the topic of Large Eddy Simulation and its 
promise for the future. 


II. EQUATIONS OF TURBULENT FLOW 
In this section, and in those that follow, we will discuss the 
steady, incompressible Navier-Stokes equations. Consider, for example, 
the continuity and x-momentum equations 


3u 3v 3w n 

37 + 37 al = 

3u , 3u , 3u 1 dp , 3 2 u , 3 2 u , 3 2 u 

3x 3y 3z p dx ^ z ^2 


dz' 


( 1 ) 

( 2 ) 


These equations contain most of the complexities associated with the 
turbulence itself, and the additional problems of compressibility, etc. 
do not serve to illuminate the discussion in any significant way. 

Suppose that the various flow parameters have time histories like 
that shown in figure 1. The velocity u can be decomposed into two 
parts as 

u(t) = u + u f (t) (3) 

where u is a steady carrier and u T (t) is a high-frequency oscillatory 
component with a zero mean. From this decomposition of the velocity we 
can define several time averages as 

■ t+T 
r t 


_ i r t+T 

u = tJ u 


(t)dt 


(4) 


and 


"'■il 


t+T 0 t+T 

[u (t) - u]dt = — J u'(t)dt = 0 


(5) 
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Similarly, 


r(t)dt ■ uv + u'v’ 


i r t+T 
fj t u(t)u 


(t)dt “ u 2 + u* 


If we substitute these composite variables into the equations of 
motion, expand, and time-average, we obtain equations like the following 


- 3u - 3u 
U "9x + V *9y + 


+ v /V U + 3fu + 

P dx \3x 2 3y 2 3z 2 / 

/Su 77 3u T v r 3uV[\ 

\ 3x 3y 3z y 


Note that we recover the original equation, in terms of time-averaged 
velocities, the barred quantities, but with the three additional terms, 
the so-called Reynolds stresses included in the rightmost parentheses. 
These additional unknowns result from time-averaging the nonlinear con- 
vective terms and have the appearance of stresses. 

When we substitute this velocity decomposition in all three momen- 
tum equations, the stress tensor for these Reynolds-averaged equations 
becomes 

(x - u f2 ) (t - uV) (x - u f w T ) 

xx xy xz 

x, . - (x v 7 u t ) (t - v f2 ) (t - v'w 1 ) (9) 

ij yx yy yz 

(x - u T w f ) (t - w'v ? ) (x - w f2 ) 

_ zx zy zz 

where the x , x , etc, are the usual viscous stresses. The nine new 
xx xy 

unknowns, the Reynolds stresses, have been introduced by the decomposi- 
tion and time-averaging of the velocities. Because of the skew symmetry 
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of the matrix only six of the nine new unknowns are independent. If we 
attempt to increase the number of equations, for example, by taking 
successive moments of the momentum equation, the number of unknowns 
continues to proliferate and this shortfall between equations and 
unknowns has been termed the "turbulence closure problem." 

The complete Navier-Stokes equations in Reyno Ids -ave raged , mass- 
averaged form, including the Reynolds stress and turbulent heat-transfer 
terms, are presented in orthogonal tensor form and cylindrical coordi- 
nates by Rubesin and Rose [7]. 

III. SIMILARITY LAWS 

Before we turn our attention to the modeling of the Reynolds 
stresses, it is useful to consider some experimental results. One of 
the most important classes of flows, and hence one for which a great 
deal of experimental data exist, is the wall-bounded shear flow. If we 
restrict our attention to a two-dimensional steady flow, and without 
complications such as wall mass transfer or separation, virtually all 
the reliable experimental data can be correlated by the Law of the Wall 
and the Law of the Wake. These similarity laws can be written in com- 
bined form, following Coles [8] , as 

u + - — = - Zn y + + - 2 sin 2 (y + C (10) 

u k J k \2 <5/ 

T 

where 


k = 0.41 
C = 5.0 
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u ■ A /p 

T W 

<$ = boundary-layer thickness 

II « wake strength 

This expression provides a good fit to the experimental velocity data over 
a wide range of pressure gradients in terms of the three parameters, , 
II, and 6. Figure 2 shows how well this correlation works: for a strong 

adverse pressure gradient, figure 2a; a zero pressure gradient, figure 2b; 
and a strong favorable pressure gradient, figure 2c. The data are taken 
from the 1968 Stanford Conference Proceedings [7]. The Law of the Wall, 
the linear portion of the correlation on the semilog plot, represents 
only about the inner 20% of the boundary layer in an adverse pressure 
gradient and about the inner 50% of the boundary layer in a favorable 
pressure gradient. The important point to note is that, despite the 
apparent chaos in the instantaneous values of the flow parameters, the 
time averages can be correlated by a judicious choice of variables; 
further, these averages are reproducible in different wind tunnels and 
with different instrumentation as long as the appropriate dimensionless 
parameters are held fixed. In addition, these correlations of experi- 
mental data provide insight into those parameters that govern the flow 
field. Such insight is critical to the construction of rational models 
of the turbulent mixing process. It should always be borne in mind that 
no improvement in turbulence modeling can take place without continued 
careful experimental studies being carried out, and without exhaustive 
analysis of the data resulting from these studies. 
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IV. TURBULENCE MODELING 

Turbulence modeling concerns itself with the generation and testing 
of closure relations describing the Reynolds stresses. Since its incep- 
tion, the goal of turbulence modeling has been the Universal model.” 

We may define this as an equation, or system of equations, capable of 
accurately describing the Reynolds stresses of any turbulent flow with- 
out any previous experimental information. Such a model does not as yet 
exist and, in fact, may never exist. At present we can confidently pre- 
dict only a relatively small class of turbulent flows or, more properly, 
a few small classes of flows. These predictable flows fall into the 
category of equilibrium flows, which in turn may be defined as flows in 
which the production and dissipation of turbulence energy, shear stresses, 
etc. are in balance. This balanced condition implies that there exists 
a one-to-one relationship between the Reynolds stresses and the mean 
flow. Physically, this means that the turbulent flow is independent of 
its history, or that it has existed for a significant time under the 
effects of mild pressure gradients and without recent or rapid changes 
in the boundary conditions, for example, blowing, bleed, slip, etc. 
Examples of such flows are boundary layers on mildly curved bodies at 
subsonic speeds, shock-free supersonic flows, far wakes, and mixing 
layers well downstream of their initiation. 

Turbulence models can be divided into three major categories: 
models suitable for integral methods, eddy-viscosity models, and 
Reynolds-stress models. We will discuss each of these categories in 
the following paragraphs. 
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A. Integral Method Models 

Many of the practical calculation schemes, including integral 
methods, are based on a subset of the Navier-Stokes equations termed the 
boundary-layer equations. For flows in which v/u is everywhere a 
small quantity, this system of equations, although substantially 
simpler than the Navier-Stokes system, contains all the relevant 
physics. In subsequent sections of the text we will consider these 
equations almost exclusively, since they simplify the discussion with- 
out significant loss of generality. The fundamental simplifying 
assumption is that variations with respect to x are much smaller than 
those with respect to y. This reduces the number of equations by one, 
for the two-dimensional case, by deleting the y-momentum equation, thus 
making the pressure distribution a parameter of the problem rather than 
part of the solution and changing the remaining system of equations to 
parabolic in x rather than elliptic. This latter change permits sub- 
stantial savings in computing costs since the resulting system can be 
solved by a marching procedure rather than by a sweeping method. The 
earliest methods developed for the computation of boundary-layer flows 
were integral methods, and such methods are still in use today. The 
principal use of these methods today is for design studies where many 
sets of calculations for the same kinds of flows are required. In this 
application where the calculation is used to interpolate between experi- 
mentally established bench-marks, integral methods are without peer. 

These methods are computationally very fast and permit the easy incorpora- 
tion of empirical information via data correlations. In discussing these 
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methods, we are called upon to discriminate between the requirements 
of engineering application and of science. In the first case, we are 
concerned with the effects of fluid flow on an object of arbitrary 
shape; in the second, we are concerned with the effects of an object of 
arbitrary sihape on a fluid flow. The choice of viewpoint determines 
the type and detail of information sought. 

In applied engineering one is concerned with such parameters as 
pressure distribution and skin friction, and any method which will pro- 
vide this information accurately and cheaply for the particular flow 
configuration at hand, is satisfactory. Given an adequate data base 
and a limited range of performance requirements, integral methods are 
quite satisfactory. 

One can obtain any of several integral methods by simply integrat- 
ing the momentum and continuity equations with respect to y. Any text 
on classical boundary-layer theory will provide examples of how this 
can be carried out. Additional relations may be obtained by taking 
moments of this equation with respect to u and its powers before 
integration (ref. [9]). 

The basis of most integral methods currently in use is Head's 
Entrainment Method [10]. In this method the momentum and continuity 
equations are integrated across the boundary layer to yield the Kannan 
Integral Equation 




A dU o 

H) - 2 -— 2 . 

u dx 
e 



(id 
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where 


5* 


•ftK) 
C (* ■ t) " 


dy 


H » — 
H 0 


w 


P U 
e e 


This provides us with one ordinary differential equation in the three 
unknowns, skin friction coefficient momentum thickness 0, and 

shape factor H. A second equation is obtained by integrating the con- 
tinuity equation across the boundary layer over the same limits, and 
combining the result with the definition of displacement thickness, 6*; 
that is , 


d ,, 5 - <S* du e 

to (6 " 6 * ) = v(6) --T— to” 


( 12 ) 


where v(5) - v e / u e is proportional to the rate at which free-stream 
fluid is entrained into the boundary layer. Using Head's notation, let 
v(S) = F and (5 - 6*) = A, we can define a new shape factor as 

H - A 
H l " 0 

where, for power-law velocity distributions, u/u^ = (y /6) n , 

H = 2H 
H l H - 1 

Green [11] notes that experimental data for Head's F(H X ) can be well 
correlated by the linear function of H 
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F = 0.025H - 0.022 


(13) 


Differentiating equation (13) with respect to x and combining 
equations (11) and (12) yield 

0 g - H(H 2 - 1) y- ^ C (H - 1)F - HC f ] (14) 

e 

Equations (11) and (14), together with a skin-friction law, for example, 
the Ludwig-Tillman Law 

C f = 0.246 exp(-1.561H)Re 0 ~° * 268 (15) 

provide the necessary three equations. The system of ordinary differ- 
ential equations is then solved by any standard integration routine, 
for example, Runge-Kutta, or Adams Moulton. For examples of these 
schemes any standard text in boundary-layer theory or the 1968 Stanford 
Conference Proceedings should be considered. The turbulence information 
is contained in equations (14) and (15) which reflect correlations of 
experimental data. 

In reference [11], Green extended this method to compressible 
flows, and in reference [12], Green and his coworkers incorporated an 
ordinary differential equation, based on the turbulence energy equation, 
to describe the entrainment rate. This latter method, the Lag-Entrainment 
Method, is an effort to incorporate flow history into the calculation. 

It is a reasonably successful procedure for some types of flows and is 
in fairly common use today. 

The advantages of methods of this type are computational simplicity 
and the ease with which the turbulence equations can be modified to 
incorporate new empirical information, that is, by simply introducing a 
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different correlation of the skin-friction coefficient. The principal 
disadvantage is a lack of flexibility in dealing with flows even slightly 
different from those for which the correlations were developed. That is 
to say, when a new flow situation is encountered new data correlations 
obtained in like flows must be used. 

Differential methods or, more properly, difference methods intro- 
duce substantially more flexibility in their description of the varia- 
tion of the mean flow, but require more detailed information about the 
turbulence. The balance of the turbulence modeling techniques discussed 
here will be those appropriate to differential mean-flow descriptions. 


B. Eddy-Viscosity Models 

For boundary-layer flows it is generally agreed that the Reynolds 
normal stresses u l2 may be ignored in comparison with the Reynolds shear 
stresses u'v’, except perhaps in the neighborhood of the separation 
point. This agreement is based on a combination of experimental data 
and general unwillingness to complicate the problem further by consid- 
ering Reynolds normal stresses. If we argue that the turbulent shear 
stress should be defined in an equivalent manner to the laminar shear 
stress, we can propose a new turbulent viscosity. Applying this argu- 
ment to the two-dimensional, incompressible boundary-layer equations, 
defined above, we can write after Boussinesq 



where e is the so-called eddy or turbulent viscosity. Note that this 


viscosity is a property of the fluid motion and not a physical property 


of the fluid itself, that is. 
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1. ZERO-EQUATION MODELS 

By zero-equation models is meant that the model requires no partial 
differential equations to describe the eddy viscosity. The current 
most commonly used algebraic eddy viscosity model is that due to Cebeci 
and Smith [13], with subsequent improvements by Cebeci [14]. This is a 
composite model, in keeping with the experimental observations embodied 
in the similarity laws noted earlier, that is, 



(19) 


In the inner region, the length scale is taken as proportional to the 
distance from the wall, with a damping term near the wall due to 
van Driest [15], and the velocity scale is taken as this distance times 
the normal gradient of velocity, that is. 
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where 



A => 0.4y^l - e ^ ^ 
is the length scale with near boundary damping, and where 

A + * 26(1 - 11.8p + )" 1 ^ 2 

+ -v dp 

P " pu 3 dx 
T 

The term p + was introduced to fit data better in pressure gradients. 
This results in an inner eddy viscosity of the form 

■i ■ [*•* (* - *~ y+/A+ )T If l f (2o> 

In the outer region the length scale is taken as proportional to the 
boundary-layer thickness and the velocity scale is taken as the so-called 
velocity defect. These definitions can be combined as 


e 

o 


f 6 

0.0168 J (u - u)dy 


( 21 ) 


In figure 3 we show a comparison of the predicted and measured 
velocity distributions in physical variables, corresponding to the same 
cases shown in similarity variables in figure 2. It can be noted that 
presentation of data in these different forms emphasizes different 
regions of the boundary layer, and can provide additional insight into 
the physics of the flow. It can be seen that the calculations based on 
the zero-equation model provide a good representation of the zero and 
favorable pressure-gradient data while missing the adverse pressure- 
gradient results rather badly. This failure appears as a significant 
underprediction of the boundary-layer thickness coupled with an 
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underprediction of the wall-shear stress. Failure in the accurate 
prediction of flows with strong adverse pressure gradients is a major 
shortcoming of these methods. As shown in figure 4 the skin friction 
as predicted with the Cebeci-Smith model is compared with the experi- 
mental data of Simpson, Chew, and Shivaprasad [16]. This flow is 
essentially a two-dimensional diffuser with a strong adverse pressure 
gradient leading to separation. The comparison is typical of algebraic 
eddy-viscosity models incorporated in a standard boundary-layer code in 
the neighborhood of separation. In these cases, the comparison is con- 
fused by the appearance of an additional complication, the separation- 
point singularity, which renders the mean-flow equations invalid at the 
separation point. This point will be discussed further in a subsequent 
section. 

From the computational point of view, models of this type are the 
simplest models that satisfy the requirements of a robust difference 
solution to the boundary-layer equations. If such models are to be used 
within a Navier-Stokes formalism, some ad hoc assumption must be made 
regarding the location of the boundary-layer edge and the definition of 
the velocity there. These parameters lack a clear definition in any 
situation in which the normal derivative of velocity does not decrease 
raonotonically as the outer flow is approached (cf. refs. [17] and [18]). 

2 . ONE-EQUATION MODELS 

As noted above, the goal of turbulence modeling is to produce the 
universal model. To redress the shortcomings of zero-equation models 
it seems plausible to require that the scaling parameters of the 
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turbulence be based on a property of the turbulence rather than on a 
property of the mean flow. In the mid 1960s Glushko [19] took a step 
in this direction by proposing that the velocity scale should be the 
square root of the kinetic energy of the turbulence, while retaining an 
algebraic length scale. The turbulence energy equation may be written 
(ref. [20]) as 


. 2 ..!, 


3q 2 /2 , - 3q 2 /2 . — j—r 3u 3 ( . q U Jl\ , 3u i 

3t- + “t 3^7 + u i u t V pu « ~T~) ‘ VU i 17 


( 22 ) 


The various terms in this equation, from left to right, are the rate 
of increase of turbulence energy, the increase of this energy due to 
the convection by the mean flow, the production of turbulence energy, 
the transport of turbulence energy by turbulent and pressure diffusion, 
and, finally, a terra that accounts for the dissipation of turbulence 
energy into heat* In attempting to solve this equation, the first three 
terms can be allowed to stand as they are while the latter two terms, 
comprising the turbulent and pressure diffusion and the dissipation, 
must be modeled, 

Rubesin [21] gives a clear presentation of the original Glushko 
model and extends it to compressible flow. In addition, he introduces 
the elliptic terms required for consistency with a full Navier-Stokes 
solution. The Glushko model for steady mean flow, including the modeled 
terms, may be written as 


- 3q 2 /2 , _ 3q 2 /2 

u — \ + V = 

3x 3y 


3 

3y 


[1 + e(Xr)] 


3q 2 /2 

3y 


+ veS 


3u . 

: 

ij 3x , 


- Cv[l + eQr)] 



(23) 
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where 


q 2 u' 2 + y' 2 + w'‘ 
2 " 2 


the turbulence kinetic energy 


v 


ij 


- A (^1 * !!i\ 

" 2 V 9x j 9 x i/ 


the dimensionless turbulent viscosity 


the mean strain rate 


and i is the turbulence length scale. To complete the model the 
relations 


where 


r _ /q 2 /2£ 


e = H(r)ar 


Cr_ 

r 


H(r ) = < 


t - 1 - 


u 


z 

<5 


i = y/6 + 0.37 
6 12.61 

1.48 - y/6 
2.42 


0 < — <0.75 
r 

o 


0.75 < — < 1.25 
r 

o 


1.25 < 


0 < -*jr < 0.23 

O 


0.23 < J < 0.57 


0.57 < |- < 1.48 


a = 0.2 


r = 110 
o 


(24) 


C - 3.93 
X = 0.4 

were suggested by Glushko. He also proposed the boundary conditions 
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y » 0 q 2 * 0 

■■ y' + 'oo 0 

This latter condition is required by the form of the equations at large 
y. The constants and functional relations were chosen to provide agree- 
ment with experimental results for flat-plate flows. 

While this model is historically interesting, as the first real 
break from equilibrium concepts in a practical calculation scheme, it 
has a major conceptual shortcoming *— that the length scale is still an 
algebraically defined quantity defined in terms of distance from the 
wall. For more geometrically complicated situations, for example, 
corner flows, backward-facing steps, etc., additional ad hoc relations 
for the length scale must be proposed. 

From the computational point of view, this model, along with the 
other multi-equation models requires an implicit calculation scheme 
because of the large range of eigenvalues (i.e., stiffness) associated 
with these equations. Given a little care, however, serious computa- 
tional difficulties can be avoided. One problem which cannot be avoided 
here is the need to specify an additional condition on the inflow 
boundary, that is, the initial turbulence energy profile. A procedure 
for obtaining such a profile (ref. [22]) is to start the calculation as 
a laminar flow with an initial energy distribution given by 

V = ~i~ (j*) exp {i f 1 - (^) ]} (25) 

where q^ 2 */2 and y*/5 are specified constants, and let the calculation 
simulate a transition process. Beckwith and Bushnell [22] conducted a 
parametric study on the effects of varying the initial energy and found 
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that at too low a value the energy simply damps with increasing downstream 
distance. This procedure is satisfactory in the absence of any experi- 
mental data. Murphy and Rubesin [18] found that if experimental velocity 
profiles are to be matched in initiating a calculation, a slightly 
inconsistent energy profile can cause the solution to have an initial 
jump from the starting condition. This behavior is shown in figure 5 
for the flat-plate data of Wieghardt (presented in ref. [5]). 

3 . TWO-EQUATION MODELS 

A much more general eddy viscosity model can be proposed by the 
introduction of a differential equation to define the length scale as 
well as the velocity scale. This is particularly important when we 
consider flows such as those over a backward-facing step or over the 
trailing edge of an airfoil. In these cases the distance from the wall 
is not well defined and, if a zero- or one-equation model is being used, 
some ad hoc assumption must be made with regard to the length scale. 
Two-equation models provide more general models by permitting the same 
length-scale equation to be used, regardless of the flow configuration. 

The principal difficulty in proposing such a model is in arriving at the 
appropriate form for the new differential equation. Since mixing length 
is a wholly artificial concept it is necessary to identify the proposed 
length scale with some average-eddy size. Taking the position that all 
turbulence properties are probably transported in an analogous fashion, 
one can model the length-scale transport equation after the turbulence 
kinetic-energy equation. Some authors have chosen not to model the 
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actual length scale, but rather a turbulence-field parameter which is 
related to it. 

One of the earliest two-equation models and the one in most general 
use today is that due to Jones and Launder [23] . The second equation 
defines a new parameter, the dissipation, which is related to the 
length scale defined previously as l » k 3 ^ 2 /e. This model can be 
written as 


3k = __i , _3_ 

3Xj “ pT ij 3xj 3Xj 


[(•*?)%]-[••*••(«?)'] 



where 



= 0.09 exp [-2. 5/(1 + 1^/50)] 
C = 1.45 

C 2 = 2 [1 - 0.3 expO-I^ 2 )] 



These constants and functional relations are obtained by matching 
experimental data in widely differing flow configurations. Note that 
we have retained the symbol k for turbulence energy in this case, 
since this model and its descendants are generally known as k-epsilon 
models. 

As might be expected, this increase in generality is not obtained 
without some cost. In particular, the computational work required by 
these methods is substantial, and in both time-marching and steady-state 
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calculations extreme care must be taken to prevent instabilities from 
propagating. These instabilities occur near surfaces where the 
turbulence-field variables are rapidly varying. The large variations 
in the scales of the problem occurring in this region is reflected in 
the computations as stiffness in the governing equations. In a recent 
paper, Viegas and Rubes in [24] overcame much of the numerical diffi- 
culty by defining wall functions to be satisfied by all the differen- 
tial equations, for compressible flow, at small y. This procedure 
apparently makes a first-order improvement in terms of the ease with 
which these equations can be solved without degrading their accuracy. 

Some ad hoc assumptions must still be made, however, in flows where an 
unambiguous distance from the wall is not available, for example, the 
edge of a backward-facing step. 

When the proper initial conditions are available for these higher- 
order methods their accuracy is, in general, comparable to that of the 
Cebeci-Smith model for flows for which the latter works well (fig. 5). 

In addition, they can be used in cases for which simpler models would 
require many changes. The cost of these methods is significant, however, 
and one should always use the simplest method consistent with the 
results sought. For a recent review of the performance of the higher- 
order methods Marvin’s discussion in reference [25] is useful. 

V. REYNOLDS STRESS MODELS 

All of the preceding results effectively ignore the question of the 
basic validity of the eddy-viscosity concept. Experimental results 
obtained over the last 20 years or so have shown that while eddy 
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viscosity is a convenient concept from the computational viewpoint, it 
has little physical basis. This should come as no surprise since the 
fundamental assumption of the arguments based on kinetic theory of the 
eddy -viscosity concept requires that, analogous to the mean free path, 
some measure of the typical eddy size be small compared with the dimen- 
sions of the container. Since the largest eddies have a dimension of 
the same order as the boundary-layer thickness, it is not surprising 
that, under some conditions, the theory fails to provide good predic- 
tions. What is surprising is that the concept works so well as long as 
we restrict our attention to flow in zero and favorable pressure 
gradients. 

A logical next step in the hierarchy of turbulence models is the 
Reynolds stress model. We will consider three models in this category 
reflecting the development of these methods over the last 15 years. 

The earliest Reynolds stress model used routinely in computation 
was that due to Bradshaw, Ferris, and Atwell [26]. In their method it 
is assumed that the Reynolds shear stress is proportional to the turbu- 
lence kinetic energy. A form of the turbulence-energy equation is then 
solved for the Reynolds shear stress directly. As with the Glushko 
model, the diffusion and dissipation terms must be modeled. This was a 
complete break from the eddy-viscosity concept, and the direct link 
between Reynolds stresses and the mean field embedded in all the models 
discussed previously. 
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The equations for the Bradshaw model may be written as 


and 



x ~ - DIFF - DISS 
xy 3y 


pq 


2 a i 


where DIFF is the turbulent diffusion term, modeled as 

1/2 


DIFF 


(¥) 


■ v Q 

9y \p 
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and DISS is the dissipation term modeled as 

DISS = (^) J/2 


and the empirical functions L and G are shown in figure 6, Note that 
in an equilibrium condition. Production * Dissipation, the dissipation 
length scale L is exactly analogous to the mixing length, and this 
equation can be used to define it. While this model is a significant 
conceptual advance over the Glushko model it suffers a similar short- 
coming in terms of generality in that the length scale is still defined 
algebraically in terms of the distance from the wall. 

A substantially more complicated model is the algebraic Reynolds 
stress model proposed by Rodi et al. [27]. This model uses the basic 
two-equation model of Jones and Launder as a starting point, but rather 
than defining an eddy viscosity, Rodi et al. chose to define the indi- 
vidual Reynolds stresses in terms of algebraic relations between the 
turbulence energy and the dissipation, or length scale. In the limit 
of equilibrium flow this method reduces to the usual Jones-Launder model 
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since the Reynolds shear-stress equation contains the same constitutive 
relation for stress as does the original model • 

In reference [28], Launder, Reece, and Rodi propose a full Reynolds 
stress model providing differential equations for the transport of all 
six of the Reynolds stresses, which must be solved simultaneously with 
the dissipation equation. These last two methods are still in the rela- 
tively early research phase and have seen only very limited application 
as of this time. As a result their usefulness and generality remain an 
open question. For the equations describing these methods the reader 
should examine the cited references. 

VI. SEPARATED FLOWS 

Separated flows have, so far, been the most difficult aerodynamic 
flows to predict accurately because 1) reliable experimental data for 
separated flow are extremely difficult to acquire, and 2) classical 
boundary- layer calculations fail at the separation point as a result of 
singularity in the equations. The rationale for considering them here 
is that these flows provide a sufficiently hard test for turbulence 
models to permit discrimination between various models. As was pointed 
out by Goldstein [29], the singularity in the boundary-layer equations 
is associated with the specification of the streamwise pressure gradient 
at the separation point. If some parameter other than pressure, for 
example, displacement thickness, skin friction, or some bounding stream 

function, is specified the boundary-layer equations are well behaved 

(■ 

despite the presence of the separation. 
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Solutions obtained in this manner, so-called inverse solutions to 
the boundary-layer equations, have been compared with solutions to the 
Navier-Stokes equations for the same distributions of skin friction and 
have been found to reproduce all the relevant physical parameters as 
long as v/u remains small, that is, as long as the flow is "slender 
The use of inverse boundary-layer methods as test beds for turbulence 
models appears to be an economical approach to the evaluation of these 
models in the vicinity of separation. Unfortunately, however, if the 
pressure is not that which was observed experimentally, no conclusions 
with regard to the specific shortcomings of the turbulence model can be 
inferred. A subterfuge was proposed by Arielli and Murphy [30], which 
permits one to specify both the additional boundary condition required 
to avoid the separation point singularity and the strearawise pressure 
distribution. This permits the specification of streamwise pressure 
distribution while retaining the economy associated with boundary-layer 
calculations. 

Several investigators, using inverse boundary-layer methods and the 
full Navier-Stokes equations, have considered the problem of turbulent 
separation and have found that, in general, current models tend to pro- 
vide too much mixing, that is, too large a value of skin friction, as 
the separation point is approached and too little mixing or too slow a 
rate of recovery of skin friction, downstream of reattachment. In an 
important aerodynamic flow situation, the shock-wave boundary-layer 
interaction, the separation point itself is reasonably well predicted, 
but since this is an inertially dominated flow, this should not be taken 
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as evidence in favor of the turbulence model. In fact, a reasonable 
prediction of the behavior of the mean flow in this case can be deter- 
mined from inviscid considerations alone (cf. ref. [31]). 

Calculations of the flow of reference [16] (the flow in a 
converging-diverging-converging channel, simulating a subsonic diffuser, 
but forcing reattachment in the final converging section) using inverse 
boundary- layer and full Navier-Stokes procedures, and specifying the 
channel configuration rather than the free-stream velocity, predict the 
separation point to be in about the correct location, but predict a 
velocity gradient in the free stream to be about twice that observed 
experimentally. These independent calculations (refs. [32] and [33]) 
are consistent with each other and with the observation that the zero- 
equation models overpredict turbulent mixing in the approach to separa- 
tion. Simpson and Collins [34] suggest that the failure of the predic- 
tions is caused by inadequate accounting for the Reynolds normal stresses 
in the neighborhood of separation. Since it is extremely difficult to 
obtain a truly two-dimensional flow in the neighborhood of a separation 
point, and equally difficult to obtain reliable measurements in the 
near-wall region, evidence to support or reject this hypothesis is 
lacking. 

Two useful references on the general topic of separated flows are 
the AGARD Conference of Separated Flow in 1976 (ref. [35]) and the 
Project Squid Summary Report on Colloquium on Flow Separation in 1979 
(ref. [36]). 
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VII. UNSTEADY AND THREE-DIMENSIONAL FLOWS 


In considering unsteady flows we must introduce a new concept, that 
of ensemble averaging. In this case we find that averaging over long 
times will introduce a confusion between actual turbulent motion and 
at least some of the high-frequency components of the mean motion. To 
avoid this confusion, we can consider the averages to take place over 
some large number of successive realizations of the same event. For 
example, if we consider a sinusoidally varying mean flow we can take a 
data reading every time we reach the same point in the cycle, and then 
take an average of these readings over many hundreds of cycles. This 
procedure is termed ensemble or phase averaging. 

Only within the last 5 years or so has the computational capability 
become available to permit routine calculation of unsteady and three- 
dimensional viscous flows. As a result, turbulence modeling for these 
flows has been almost exclusively an adaptation of the two-dimensional 
steady models discussed above. As in the steady case, predictions of 
unsteady flows in favorable and mildly adverse pressure gradients can 
be made with reasonable accuracy as long as the dimensionless frequency is 
not too high. At higher frequencies the existing experimental data are 
inadequate to permit a firm judgment to be made. References [37] and 
[38] provide an introduction to the literature of unsteady turbulent 
flows. 

For truly three-dimensional flow (e.g., ref. [39]) too few calcula- 
tions have been made to permit any general conclusions to be drawn. An 
example of some of these calculations is the work of Kussoy, Viegas, 
and Horstmarm [40]. 
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VIII. LARGE EDDY SIMULATION 

Within the last 5 years, investigators with access to very large 
computers have been obtaining solutions to the three-dimensional, 
unsteady Navier-Stokes equation on a relatively fine mesh. These solu- 
tions, for low Reynolds number turbulent flows, have sufficient spatial 
and temporal resolution to capture at least some of the turbulent eddy 
structure. Moin et al. [41], as reported by Rubesin [42], have obtained 
solutions of a channel flow at a Reynolds number, based on channel half 
width, of 13,800. These conditions correspond to the experiment of 
Hussain and Reynolds [43]. In figure 7 (taken from ref. [42]), we show 
two velocity profiles corresponding roughly to the maxima of the posi- 
tive and negative excursions in the instantaneous velocity, together 
with the long time average of the resolved velocities. As can be seen, 
the appearance of these profiles gives a quite plausible representation 
of a turbulent flow. In figure 8 the same long time averaged profile 
plotted on the usual turbulent similarity coordinates are shown. 

Despite the relative crudeness of the calculation, and the various 
assumptions regarding initial conditions and subgrid modeling, these 
results provide substantial hope for the ability to compute turbulence 
in the long term, and in the nearer term such calculations may provide 
substantial insight into the modeling of parameters which presently 
cannot be measured. 

IX. SOME CONCLUDING REMARKS 

In the foregoing we have considered only relatively simple flows 
of the type used to evaluate turbulence models. In a situation in which 
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the reader is called upon to make use of the information contained here, 
the flow situation will almost certainly be much more complicated. Since 
one cannot provide the appropriate information for every application, it 
seems more appropriate to provide the reader with an appreciation of 
where we are with respect to turbulence modeling, and of how we got 
there. With this information and an acquaintance with the material in 
the works referenced here, it is hoped that the reader will be in a 
better position to make rational judgments when confronted with more 
complex problems. Some additional references which may be found useful 
are references [44]-[46]. 

We have not considered any of the difficulties associated with 
such complicating factors as heat transfer, mass transfer, and/or chemi- 
cal reactions since the problem of turbulence is "sufficient unto the 
day." In addition to the references cited here the author has frequently 
found the following journals to be helpful: the several Transactions of 

the ASME, the International Journal of Heat and Mass Transfer, and the 
various publications of the AIAA. Workers in specialized fields other 
than aerodynamics will undoubtedly require and be familiar with different 
sources. 
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LIST OF SYMBOLS 


F 

G 

H 

H i 

k 

L 

l 

P 

+ 

P 

q 2 /2 

Re 

T 


t 


u,v,w 


x,y,z 

+ 

y 


6 


6 * 

A 


function used in defining near-wall damping in Cebeci-Smith eddy 

viscosity relation (ref. [13]) 

skin friction coefficient, x / (l/2)p u 2 

w e o 

function defined in equation (3) 

diffusion function in Bradshaw et al. model (ref. [26]) 
boundary- layer form factor, S*/0 

Head’s boundary layer form factor, A/0 (ref. [10]) 
turbulence kinetic energy in Jones-Launder model (ref. [23]) 
dissipation length scale, model of Bradshaw et al. (ref. [26]) 
mixing length 
pressure 

dimensionless pressure gradient parameter in A + 
kinetic energy of turbulence 
Reynolds number 

a time period, long relative to the turbulent motion period, but 

short compared with mean motion period 

time variable 

velocities 

shear velocity 

spatial coordinates 

a wall unit Reynolds number 

boundary-layer thickness 


displacement thickness = J 

o 

Head’s parameter =6-6* (ref. [10]) 


1 - (Au/p e U Q )dy 
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FIGURE CAPTIONS 


FIG. 1. History of instantaneous velocity, a) Hot-wire signal. 

b) Reynolds decomposition. 

FIG. 2. Comparison of experimentally determined velocity; profile with 
the Law of the Wall-Law of the Wake, a) Strong adverse pressure gradient 
(data of Ludwig-Tillman) ; ID 1200* x =* 2.782 m. b) Zero pressure gra- 
dient (data of Wieghardt); ID 1400, x * 3.487 m. c) Strong favorable 
pressure gradient (data of Herring and Norbury) ; ID 2800, x ■ 2.0 ft. 

FIG. 3. Comparison of experimental and predicted velocity profiles, 
a) Strong adverse pressure gradient, b) Zero pressure gradient. 

c) Strong favorable pressure gradient. 

FIG. 4. Comparison of experimentally determined skin friction coeffi- 
cient with the prediction of an algebraic eddy viscosity model. 

FIG. 5. Comparison of experiment skin friction of a flat plate with 
predictions using three eddy viscosity closure models. 

FIG. 6. Length scale and dissipation functions for the model of 
Bradshaw et al. (ref. [26]). 

FIG. 7. Instantaneous velocity profiles and Reynolds averaged velocity 
profiles from large eddy simulation calculation. 

FIG. 8. Reynolds averaged velocity profile of Fig. 7 compared with the 
experimental data of Hussain and Reynolds (ref. [43]). 
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